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Abstract 

In this paper we continue the investigation of the lattice gas model. The main 

improvement is that we use two strengths for bonds: one between like particles 

and another between unlike particle to implement the isospin-dependence of 

nuclear force. The main effect is the elimination of unphysical clusters, like the 

dineutron or diproton. It is therefore a better description of nuclear system. 

Equation of state in mean field theory is obtained for nuclear matter as well 

as for N ^ Z systems. Through numerical and analytical calculation we show 

that the new model maintains all the important features of the older model. 

We study the effect of the Coulomb interaction on multifragmentation of a 

compound system of A=86, Z=40 and also for A=197, Z=79. For the first 

case the Coulomb interaction has small effect. For the latter case the effect is 

much more pronounced but typical signatures of the lattice gas model such 

as a minimum (maximum) in the value of r (S%) are still obtained but at a 

much lower temperature. 
25.70.Pq, 24.10.Pa 
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I. INTRODUCTION 



Over the past a few years we have been developing a lattice gas model for the study of 
nuclear multifragmentation[l-5]. In this model n nucleons are placed in N cubes and they 
interact with nearest neighbour interactions. In most of this work the interaction between 
neutron-neutron, proton-proton and neutron-proton was taken to be identical although we 
have sometimes used [5] a more complicated model in which interactions between like parti- 
cles are different from those between unlike particles. The purpose of this paper is to more 
fully expose this improved model and to examine its relationship with our earlier and simpler 
model. Among other things we will also show that the important conclusions reached with 
our simpler model go through in our improved version. 

The lattice gas model has the attractive feature that it can provide an equation of state 
but it can also provide the cluster distribution. Most of the observables are calculated 
using Monte-Carlo simulations. We used Metropolis algorithm in our simulations. The n 
particles are distributed in N boxes according lattice gas Hamiltonian and their momenta are 
generated from Maxwell-Boltzmann distribution at a prescribed temperature. Calculation 
of clusters is straightforward. Two nucleons in neighbouring cells are part of the same cluster 
if the relative kinetic energy is less than the strength of the attractive bond: p 2 r /2[i + e < 
0. Here p r is the relative momentum, (j, the reduced mass and e is negative (attractive 
interaction). This prescription is sufficient to calculate cluster distribution. The motivation 
for introducing two kinds of bond is now obvious. If the neutron-neutron bond or the 
proton-proton bond is attractive then each numerical simulation can generate dineutrons 
and diprotons which in reality do not exist in nature as composites. One can get rid of these 
unphysical clusters by simply making the neutron-neutron and proton-proton bonds zero or 
repulsive. 

For completeness we mention some more features of the lattice gas model that were 
studied in [1,2]. One finds that at a certain temperature the distribution of composites is 
a power law: Y(Z) oc Z~ T where Y(Z) is the number (averaged over many simulations) of 
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composites with Z protons. As is a common practice, we will extract a value of r even when 
the distribution has significantly deviated from a power law[6]. This value is extracted by 
using 

Y.fZY(z) _ zI°zz-t 
El°Y(z) El°z- 

We also calculate the second moment S 2 = J2' A 2 Y(A)/n where in the sum the largest 
cluster is excluded. The usefulness of S 2 was emphasized by Campi[7]. Since the lattice gas 
model has a Hamiltonian, its average energy at any temperature can be calculated. The 
excitation energy and specific heat C v per particle can be calculated. We will find this useful 
too. 

II. THE EQUATION OF STATE 

We have N boxes in the lattice in which we have to put n n neutrons and n p protons; 
n p + n n = n < N; n/N = V /V = p/ p where V , p are normal nuclear volume and density 
respectively. In principle, we can have three kinds of bonds: e nn , e pp and e pn . In the most 
general case where these interactions can take any arbitrary values a rich assortment of 
phenomena is predicted. The grand canonical partition function of this general lattice gas 
model can be mapped on to a spin 1 Ising type model in the presence of a magnetic and 
quadrupole field. These have been studied in detail in the past [8]. For the nuclear case the 
values of the interactions are quite restricted and the richness of phenomena disappears. First 
of all we have to set e nn and e pp to be either zero or repulsive so that one avoids producing 
unphysical dineutron or diproton bound clusters. Charge independence of nuclear forces 
suggests that we put e nn = e pp . From now on we will write e pp for both e pp and e nn . In our 
past work [5] and also in other modelings [9] of nuclear collisions using classical mechanics a 
slightly repulsive e pp was used. To avoid proliferation of parameters we will set this bond to 
zero in this work. The binding energy of nuclear matter fixes the value of e pn at -5.33 MeV. 

Throughout this work 7 stands for the number of nearest neighbours. In 3 dimensions 
one has 7 = 6. We use the Bragg- Williams mean field theory using the canonical ensemble. 
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There are N boxes and n p protons and n n neutrons. Let one of the boxes be occupied 
by a proton. Then, in the Bragg- Williams approximation, among its nearest neighbours, 
on the average r yn p /N will be occupied by protons and 7n n /iV by neutrons. The number 
of n — p bonds will be 7n p n n /iV, the number of p — p bonds will be (l/2) r yn p n p /N where 
the factor of 1/2 remedies the double counting for proton-proton bonds. Similarly starting 
with a box occupied by a neutron we come up with the same number of neutron-proton 
bonds and the number of neutron-neutron bonds is determined to be (l/2) r yn n n n /N. Thus 
the interaction energy when there are n p protons, n n neutrons placed in N boxes is E = 
r y(e pn n p n n + e nn {v? n + n 2 )/2)/N and the partition function is 

z ^ n ^ = W^hw. eM ~ m (2 ' 1) 

We now find pressure P from the equation P = kT[d\n Z/dVjr and V = a 3 N where 
a 3 = 1/po is the volume of each box. Using Stirling's formula one arrives at 

P = Po kT\nj^— + p Ql e pn {n p /N){n n /N) + Po ie nn {{n n /N) 2 + (n p /N) 2 )/2 (2.2) 

Introduce an asymmetry parameter r\ = (n n — n p ) / (n n + n p ) which takes value 1 for neutron 
matter, for nuclear matter and -1 for proton matter. We can then write 

Determine the critical point from dP/dp = d 2 P/dp 2 = 0. This gives the critical density 
p c = .5po an d the critical temperature to be — / 4)[(e pn + e pp ) / 2 + (1 /2)r] 2 (e nn — e pn )}. In this 
approximation with e pn attractive and e nn = the critical temperature for nuclear matter 
(which has rj—0) would be highest at — (7/4)(e pn /2) and would fall off quadratically with rj 
to at neutron or proton matter. 

The Bragg- Williams approximation is the simplest mean-field approximation. An im- 
proved treatment using Bethe-Peierls approximation is worked out in the appendix. The 
mean field calculation shown in this section and the appendix is merely to form a rough 
idea about the nature of phase transition. In practical calculations we need to obtain the 



yields of the composites at a given temperature. Mean field theories do not provide these 
and we need to do event by event calculation which can be obtained through Monte-Carlo 
samplings. 

III. MONTE-CARLO RESULTS 

As far as we know exact results with two kinds of bonds are not available. For one kind of 
bond one can often interpret essentially exact although numerical results from well-studied 
spin 1/2 Ising model for use in the lattice gas model. In the absence of such exact results 
our only recourse is to compare numerical results obtained with two kinds of bonds, i.e., 
e P n = —5.33 MeV, e nn = with those obtained with one kind of bond that we have used 
before, i.e., e pn = e nn = —5.33 MeV. We use here N = 7 3 , which is a number appropriate 
for finite systems that we will investigate. As mentioned in the introduction the calculation 
proceeds by first putting the required number of protons and neutrons in the N boxes using 
a standard Metropolis algorithm. Nucleons are then assigned momenta from Monte-Carlo 
sampling of a Boltzmann distribution at the given temperature. The energy of the event 
can now be calculated. Clusters are then determined as explained in the introduction. 
The results shown in Figs. 1 and 2 are obtained by averaging over 1000 events for selected 
temperatures. For two assumed freeze-out densities we calculate the specific heat, the second 
moment S2 and the deduced values of r. The unit chosen in the graph for temperature is 
T c = 1.1275|e pn | which is the T c for an infinitely large lattice with one kind of bond. We find 
that the peaking of C V ,S2 and the minimum of r happen at a slightly lower temperature 
(about ten percent) with two kinds of bond as compared to when the same bond is used 
for all the particles. For example for the minimum of r to appear at the same temperature 
the value of e pp = e pn has to be set at about 10% lower value than the value of e pn when 
e pp is set to zero. Qualitatively and even semi-quantitatively the results in the two models 
look similar when this renormalisation of the strength is done. However, the peaking of C v 
is more pronounced in the two bonds model. 
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All models which employ freeze-out densities assume that the freeze-out density is less 
than .5p - If the freeze-out density is less than .5p then in the lattice gas model a peak in 
the C v will signify the crossing of the co-existence curve and a first order phase transition. 
The value of specific heat can be deduced from the caloric curve [10] but locating the peak is 
very difficult in experiment. In a recent paper we have suggested [11] that since the peaking 
of C v is accompanied by a minimum in r and a maximum in S2, the appearance of the last 
two could be taken as a signal of the phase transition. The appearance of the maxima in 
S2 and of the minimum in r in close vicinity of the maximum in C v happens in both the 
versions of the lattice gas model. 

IV. A STUDY ON THE EFFECTS OF THE COULOMB FORCE 

Here we follow the methods employed in ref [3]. At that time we studied the influence of 
the Coulomb force on fragmentation of a system which had 85 nucleons and found the effects 
to be small. For a much larger system (Au on Au: central collision so that the compound 
system has A 394) we found the effect to be very large. One of the rather unavoidable 
features of the lattice gas model is the appearance of a minimum in the extracted value of r 
as a function of temperature. This feature disappeared for the very large system of A = 394 
because of the Coulomb interaction. 

For completeness a short description of a similar calculation but done for two kinds of 
bonds will be given here. In addition to lattice gas calculations, we do molecular dynamics 
calculations whose purpose is two-fold. One is to check if the predictions of a lattice gas 
model can resemble those of a molecular dynamics calculation provided the initial conditions 
are the same and the forces are chosen to be such that they resemble implied forces of the 
lattice gas model. For this we place the n p protons and n n neutrons in the N boxes using 
as usual Metropolis algorithm. Next we assign the momenta from Monte-Carlo sampling 
of a Maxwell-Boltzmann distribution. Once this is done the lattice gas model immediately 
gives the cluster distribution using the rule that two nucleons are part of the same cluster 
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if pJ;/2/! + e < 0. To calculate clusters using molecular dynamics we propagate the particles 
from this initial configuration for a long time under the influence of the chosen force (we will 
give the force parameters shortly). At asymptotic times the clusters are easily recognised 
(a detailed discussion of cluster recognition which requires shorter computer times can be 
found in ref 12). The cluster distribution in the two models can now be compared. Fig. 3 
shows the two prescriptions give nearly the same answer. 

We now come to the second and more important purpose of the molecular dynamics 
calculation. We now add the Coulomb interaction to the nuclear part. The initialisation of 
putting the nucleons in iV boxes is done again but now with the inclusion of Coulomb forces. 
We then do a molecular dynamics propagation including the Coulomb force. The clusters 
can again be calculated and compared with the cases where the Coulomb force was ignored. 

We now give the force parameters for molecular dynamics propagation. The neutron- 
proton potential was taken to be v pn (r) = A[B(r /r) p — (r /r) q ] exp([l/(r/r — a)]) for 
r/r-Q < a and v pn (r) = for r/r > a. Here r = 1.842 fm is the distance between the 
centers of two adjacent cubes. We have chosen p — 2, q — 1, a — 1.3, B = .924 and A = 1966 
MeV. With these parameters the potential is minimum at r with the value -5.33 MeV, is 
zero when the nucleons are more than 1.3r apart and becomes strongly repulsive when r 
is significantly less than r®. We now turn to the nuclear part of like particle interactions. 
Although we take e pp = in lattice gas calculations the fact that we do not put two like 
particles in the same cube would suggest that there is short range repulsion between them. 
We have taken the nuclear force between two like particles to be the same expression as 
above plus 5.33 MeV upto r = 1.842 and zero afterwards: v pp (r) = v pn (r) — v pn (r ) for 
r < r and afterwards. This means there is a repulsive core which goes to zero at r and 
is zero afterwards. 

The results shown in Figs. 3 and 4 can be summarised as follows. Fig 3. first of all 
shows that if there is no Coulomb interaction then lattice gas model results are quite close 
to that of molecular dynamics simulation provided in the latter one starts from the same 
initial condition and uses a force suitably chosen. Fig 3 also shows that in the case of 
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A = 85, Z = 40 the Coulomb force does not have a large effect. The minimum in r and the 
maximum in £2 are shifted to slightly lower temperature. The effect for A = 197, Z = 79 is 
much bigger. The minimum in r and the maximum in £2 are shifted from 4.8 MeV (lattice 
gas without Coulomb) to about 2.4 MeV. Our previous calculation showed that there is no 
minimum in r for A = 394, Z = f 58. So somewhere between these two limits the minimum 
will vanish. 

V. DISCUSSION 

The two bond model is a natural progression of the simpler lattice gas model. In this 
paper we have done calculations with the two bond model. Although in detail the two 
models differ the major characteristics of the well studied simple model remain unchanged. 
The lattice gas model remains a quick tool to calculate experimental data. 

When the Coulomb force is very strong the lattice gas model can not be relied upon. Figs. 
3 and 4 give some indication of the reliability of the model in the presence of a Coulomb 
force. Calculations above indicate that a viable (although much more time consuming) 
prescription might be: obtain the initial conditions as in a lattice gas model. Put n nucleons 
in N boxes by Metropolis sampling where one includes in addition to lattice gas Hamiltonian 
the Coulomb force. Obtain momenta of each nucleon from a Monte-Carlo sampling of a 
Maxwell-Boltzmann distribution. Then propagate by molecular dynamics to obtain cluster 
distributions. Techniques developed in ref. 12 might be useful so that one does not need to 
run molecular dynamics till asymptotic times. One attractive feature of this hybrid model 
is that the Coulomb force is operative even during the formation of clusters as opposed to 
other models where the Coulomb force only adds repulsion between the composites already 
formed. 
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APPENDIX A: 



We follow the method of reference 1. We break up the lattices into N/(j + 1) blocks, 
each of which contains 1 central box and 7 nearest neighbours to it. We refer to Fig A.l 
where for simplicity a two-dimensional lattice is shown. The interactions within each block 
are taken into account exactly while the interactions between different blocks are treated in 
an approximate fashion. The grand partition function can be written as the product of the 
grand partition functions of the N/ (7 + 1) blocks: 

N 

Z gr = z gr (block!) z gr (block2) z gr (block ) (Al) 

7 + 1 

We want to write down the grand partition function of the block denoted by 1,2,3,7 an d 
5. In the general case there will be two absolute fugacities \ p and A n , the first referring to 
protons and the second to neutrons. The grand partition function for a block can be written 

as 

z gr = A + B + C (A2) 

where 

A = (1 + e^-^ + e^-^y (A3) 
B = e Ap (l + e Ap ~ /3fp ~ /3epp + e ^n-/3e n -Pe pn y (A4) 

C = e A,l (l + ^v-^v-^v^ _|_ e ^n-Pe n -fSe pp y (A5) 

Here A=contribution to the partition function when the central (innermost) site is 
empty, 5=contribution to the partition function when the central site has a proton and 
C=contribution to the partition function when the central site has a neutron. Equation (4) 
takes the place of eq.(3.9) in ref [1]. The two constants e p and e„ are the average interaction 
energy with the adjacent block when a proton (neutron) occupies a peripheral site. 
The probability that the central site is occupied by a proton is n p /N. Thus we have 
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(A6) 

But since no particular site is favoured over another one, the average occupation occupation 
of one of the peripheral sites must also be n p /N. This gives 

n p _ E + F + G 



TV Zg T 



(A7) 



where 

E = e Ap -^(l + e Xp -^ + e^-^y- 1 (A8) 

p — g^PgAp-/3ep-/3e pp ^ _|_ e \- /3e p - /3e pp _|_ e \„-(3e„-/3e pn ^-y-l (A9) 

Q — e ^n e Xp-/3e p -/3e pn ^ _|_ e \ p -l3t p ~-f3e pn _|_ e X„-fBe„~/3e pp y-l (AlO) 



Similarly two equations can be written for n n /N. 

n n _ C 

N Zg r 



nn _ I + J + K 

N Zgj, 



(AH) 



(A12) 



where /, J, K can be written down from expressions E, F, G by interchanging protons with 
neutrons. 

Equations (8),(9),(13)and (14) determine the four constants X p , X n ,e p ,e n . 
For n p = n n , the calculations simplify. Now we have A = A n = \ p and e = e n = e p . Then 
B = E + F + G (eqs(8) and (9)) leads to 

e A (l + Qe^y- 1 = e x -^{l + 2e A - /3e ") 7 " 1 + e 2A (l + Qe^y^e^Q (A13) 

where we have defined Q = e _/3pp +e _/3pn . Dividing both sides of eqn (15) by e x (l+Qe x ~ /3s ) 1 ~ 1 
we obtain 

l = e-*'( + T\ „- ] (A14) 
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Rewrite eq.(8) as 2N/n = z gr /B where n = n p + n n to obtain 

2N „ . _ A ( 1 + 2e A -^\ 7 



= 2 + e i ^ i a- ( Al5 ) 

Using eq(15) the above relation leads to 

e A = n e ^-7/(7-D (A16) 
2(iV - n) 

Define a; = e -^/(T-i). Going back to eq(17) one can now derive a simple solution for x: 

1 



X= 2 



N-2n (N^2n\\ ^ 
N-n \\N-nJ ^N-n 



(A17) 



Values of e and e A from the definition of x and eq.(18). 

Let us now go back to the partition function for the lattice as given by eq.(3). where 
it is written as a product of the partition functions of the N/(j + 1) blocks. If we simply 
z gr for each little block as calculated above we will count twice the interaction between 
neighbouring sites in different blocks. For example, the binding energy between 1 and 6 
(Fig. A.l) is included in z 9r (block 1) and it is included again in z 9r (block 2). We note that 
on the average there are n/N particles at each site and each block has 7 peripheral sites. 
Thus when we evaluate the partition function for the lattice, the partition function for each 
block should be corrected by the multiplicative factor 

correction = e ^ n/N (A18) 

We can now use PV = kT In Zgr, V = N/ p , and In Z gr = N/ (7 + 1) In z gr to obtain 



1 

7T1 



P = pq/cT— — j- In z gr (A19) 



where z gr includes the correction factor. 



In Fig. A2 we have drawn P — V diagrams for nuclear matter for two cases: (1) e pn = 
£ P p = —5.33 MeV and (2)e pn = —5.33 MeV and e pp = MeV. In Bethe-Peierls approximation 
the T c in the former case appears to be between 6 and 8 MeV and in the second case between 
4 and 6 MeV. 
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FIGURES 

FIG. 1. A comparison of the calculated values of r, C v and the second moment 52 in the one 
type of bond model (top panel) and two types of bonds model (bottom panel). Here as elsewhere 
T c = 1.1275|e p „| w 6 MeV. The compound system has A=103,Z=45. Notice the maxima and the 
minimum shift to lower temperature when e pp is set to zero. Also C v is more sharply peaked. 

FIG. 2. The same as above except that a higher freeze-out density is used. The number of 
lattice sites is still 7 3 . Here A=171, Z=70. 

FIG. 3. The top part compares the r values extracted from a lattice gas calculation (dotted 
curve) with those extracted from a molecular dynamics calculation which had no Coulomb (dashed 
curve) and one which had Coulomb included in the molecular dymamics calculation (solid curve). 
Molecular dynamics without Coulomb gives results very similar to those of the lattice gas model. 
Here the effect of the Coulomb is small. The number of protons was 45. The lower part compares 
the second moments. 

FIG. 4. Effect of Coulomb on r and the second moment for a much larger system; A=197, 
Z=79. The minimum in r and the maximum in S% shift from about 5 MeV to 2.4 MeV because of 
Coulomb effect. At some larger Coulomb field the minimum in r will finally disappear. 

FIG. Al. A square lattice is divided into blocks to illustrate the Bethe-Peierls approxi- 
mation. See text for details. 

FIG. A2. p-V diagram in the Bethe-peierls approximation with same bond strength -5.33 
MeV (top panel) between all particles and for the case where we distinguish between like 
particle interaction and unlike particle interaction (bottom panel). 
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